#/opt/software/R/3.1.1/bin/R

# make_plots.R
# Caleb Matthew Radens
# cradens@mail.med.upenn.edu
# 2016_1_6

# This script plots the top results from wrapper.R (identifies high hypothesis 4
#  genes from summary_table.txt, and then runs coloc.abf on them, then plots them.)

# I tested this script using the PMACS R module 3.1.1
system("echo ===========================",wait=FALSE)
system("echo inside make_plots.R",wait=FALSE)
R_ver <- substr(version$version.string,1,15) # Get R version
system(paste("echo",R_ver),wait=FALSE)
system("echo ===========================",wait=FALSE)

if (substr(R_ver,11,15) != "3.1.1"){
  stop("Please load R module 3.1.1 before initiating this script")
}

# Choose a path to load R packages from:
lp <- "/project/chrbrolab/analysis/cradens/bin/r_libs/r_module_3_1_1"

# Add your library path to the current session of R's library path variable
.libPaths(lp)

source("import.R")
source("coloc_analysis.R")
source("plot_coloc.R")

system("echo wrapper package and script dependencies loaded and checked",wait=FALSE)

# Directories in chrbrolab:
GWAS_directory <- "/project/chrbrolab/GWAS/"
eQTL_directory <- "/project/chrbrolab/CLI/snptest_output/"

if (length(list.files(pattern = "summary_table.txt"))==0){
  stop("Please make sure summary_table.txt is in the same directory as make_plots.R or modify this script accordingly.")
}

# Import summary table generated by wrapper.R
table <- read.table("summary_table.txt",header = FALSE, stringsAsFactors = FALSE,colClasses = 
                      c("character",
                        "integer",
                        rep("double",5)))

#Extract genes and traits for the rownames and add them to the table as their own columns
#(Also, remove the TES from the rownames)

names(table) <- c("gene_trait_tes", "nsnps", "hyp0", "hyp1", "hyp2", "hyp3", "hyp4")
names<-table[,1]
names<-matrix(unlist(strsplit(names,split="_")),ncol=3,byrow=TRUE)
gene_full <- names[,1]
# Note: a period is a regular expression. 
#   - two backslashes escapes the period's meaning.
gene <- matrix(unlist(strsplit(gene_full,split="\\.")),ncol=2,byrow=TRUE)[,1]
traits <- names[,2]
names<-paste(gene_full,traits,sep = "_")
table$gene_trait <- names
table$trait <- traits
table$gene <- gene

# Remove all rows with hyp4 less than some cutoff probability, sort table by hyp4
cutoff_probability <- 25
table$hyp0 <- round(table$hyp0*100,1)
table$hyp1 <- round(table$hyp1*100,1)
table$hyp2 <- round(table$hyp2*100,1)
table$hyp3 <- round(table$hyp3*100,1)
table$hyp4 <- round(table$hyp4*100,1)
table_copy <- table
table <- table[which(table$hyp4>cutoff_probability),]
table <- table[with(table, order(-hyp4)), ]
table_copy <- table_copy[with(table_copy, order(-hyp4)), ]

my_genes<-table$gene

# The following is a list of genes identified by giambartolomei in her coloc paper
# as well as the genes pointed out by the GLGC paper

# giamb_genes <- c("SDC1","TGOLN2", "INHBB","CMTM6","C6orf106","UBXN2B","VLDLR","VIM",
#                  "CYP26A1","CUX2","OGFOD1","HP", "HPR", "PPARA")
# giamb_genes <- c(giamb_genes,"CELSR2","PSRC1","SORT1", "PSMA5", "TGOLN2","SLC39A8","PPP1R3B",
#                  "TTC39B","FADS1","MMAB","CUX2","ALDH1A2","LIPC",
#                  "LIPG","LILRB2","LILRA3","ANGPTL4","PLTP")
# giamb_genes <- c(giamb_genes,"CELSR2","PSRC1","SORT1","PSMA5","ANXA9","TMEM57","INHBB","UBXN2B",
#                  "PPP1R3B","VLDLR","ST3GAL4","SPTY2D1","FADS1",
#                  "CUX2","NYNRIN","HP","HPR","KPNB1","SPTLC3")
# giamb_genes <- c(giamb_genes,"GCKR","C2orf16","CYP26A1","FADS1","LIPC","VKORC1",
#                  "HP","HPR","PLTP")
# GLGC_genes <- c("PIGV","RRNAD1","C1orf220","CPS1","ATG7","SETD2","RBM5","STAB1","GSK3B","C4orf52",
#                 "FAM13A","ADH5","RSPO3","DAGLB","SNX13","IKZF1","ABP1","MARCH8","OR4C46","PCNXL3",
#                 "MOGAT2","ZBTB42","FPR3","ANXA9","EHBP1","INSIG2","LOC84931","FN1","CMTM6","DNAJC13",
#                 "CSNK1G3","MIR148A","SOX17","BRCA2","APOH","SPTLC3","SNX5","MTMR3","ASAP3","ABCB11",
#                 "FAM117B","UGT1A1","PXK","KCNK17","HBS1L","C7orf50","VLDLR","VIM","PHLDB1","PHC1",
#                 "DLG4","TOM1","PPARA","DOK7","VEGFA","MET","AKR1C4","PDXDC1","MPP3","INSR","PEPD")
# genes<-unique(c(giamb_genes,GLGC_genes,"CSNK1G3","GPIHBP1","RBM6","RHCE","ACP2"))

# all_genes <- get_gene_names("/project/chrbrolab/CLI/snptest_output")

traits <- c("LDL","HDL","TC","TG")

results <- batch_run_coloc_abf(eQTL_directory, GWAS_directory, genes=my_genes,traits=traits,summary_only = FALSE,significant_only = FALSE)

print("coloc analysis finished for the following gene-by-traits:")
print(names(results))

# Range indicates flanking distance from gene
# X_axis_ticks indicates how many ticks should be on the final plot
# The plots are all outputted into the current working directory
plot_coloc_results(results,Range=2e5,X_axis_ticks=13)